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Abstract 

We discuss simplified models for photo-meson production in cosmic accelerators, 
such as Active Galactic Nuclei and Gamma-Ray Bursts. Our self-consistent models 
are directly based on the underlying physics used in the SOPHIA software, and can 
be easily adapted if new data are included. They allow for the efficient computation 
of neutrino and photon spectra (from ir° decays), as a major requirement of modern 
time-dependent simulations of the astrophysical sources and parameter studies. In 
addition, the secondaries (pions and muons) are explicitely generated, a necessity if 
cooling processes are to be included. For the neutrino production, we include the 
helicity dependence of the muon decays which in fact leads to larger corrections than 
the details of the interaction model. The separate computation of the ir°, tt + , and tt~ 
fluxes allows, for instance, for flavor ratio predictions of the neutrinos at the source, 
which are a requirement of many tests of neutrino properties using astrophysical sources. 
We confirm that for charged pion generation, the often used production by the A(1232)- 
resonance is typically not the dominant process in Active Galactic Nuclei and Gamma- 
Ray Bursts, and we show, for arbitrary input spectra, that the number of neutrinos 
are underestimated by at least a factor of two if they are obtained from the neutral 
to charged pion ratio. We compare our results for several levels of simplification using 
isotropic synchrotron and thermal spectra, and we demonstrate that they are sufficiently 
close to the SOPHIA software. 
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1. Introduction 



Photohadronic interactions in high-energy astrophysical accelerators are, besides proton-proton 
interactions, the key ingredient of hadronic source models. The smoking gun signature of these in- 
teractions may be the neutrino production from charged pions, which could be detected in neutrino 



telesc opes (jANTARES Collaboration 



1999 



Karleetal 



20051 ). such as IceCube; see, e.g., iRachen Sz Meszarosl (|1998l ) for the general theory of the astro 



200d;lTzamariasl2003l : lNEMO Collaboration 



physical sources. For instance, in GRBs, photoh adronic interactions are e xpected to lead to a 
significant flux of neutrinos in the fireball scenario (jWaxman fe Bahcalllll997l ). On the other hand, 
in AGN models, the neutral pions produced in these interactions may describe the second hump 
in the observed photon spectrum, depending on the dominance of synchrotron or inverse Compton 
cooling of the electrons. The protons in these models are typically assumed to be accelerated in the 
relativistic outflow together with electron and positrons by Fermi shock acceleration. The target 
photon field is typically assumed to be the synchrotron photon field of the co-accelerated electrons 
and positrons. Also thermal photons from broad line regions, the accretion disc, and the CMB 
may serve as target photon field. The latter case, relevant for the cosmogenic neutrino flux, is not 
considered in detail. 

T he basic ideas of complete hadronic models for AGN have been described already in previous 



works ( Mannheim 



1993 



Miicke fc Protheroe 2001 



Aharonian 



2002), as we ll as leptonic models 



have been discussed by iMaraschi et al.l (119921 ): iDermer fc Schlickeiserl (119931 ). In the first place, 
these models have been used as static models to describe steady-state spectral energy distribu- 
tions. But with today's generation of gamma-ray telescopes, a detailed analysis of the dynamics of 
very high energetic sourc es is possible, and time-dependent model i ng; is inevitable. For the case of 



leptonic models, see, e.g.. iBottcher &; Chiangj (120021 ) : iRiiger et al.l (|2010l ). A necessary prerequisite 



for a time-dependent hadronic modeling is the efficient computation of the photohadronic interac- 
tions. The on-line calculation via Monte C arlo simulations is not feasib le, therefore a parametric 
description is the most viable way; see, e.q.. lKelner fc Aharonian! (|2008l ). 



The prediction of neutrino fluxes in many source models relies on the photohadronic neutrino 
production. In this case, astrophysical neutrinos are normally assumed to originate from pion de- 
cays, with a flavor ratio at the source of (/ e , / M , f T ) ~ (1/3, 2/3, 0) arising fro m the decays of both 
prim a ry pions and secondary muons. However, it was pointed out in Ref. (jRachen fc Meszaros 
199a ; iKashti fc Waxmanl 120051 ) that such sources may become opaque to muons at higher ener- 
gies, in which case the flavor ratio at the source changes to (/ e , / M , f T ) ~ (0, 1, 0). Therefore, 
one can expect a smooth transition from one type of source to the other as a fun ction of the 
neutrino energy (|KachelrieB k, Tomasl 12000 ; iLipari et al.l 120071 ; iKachelriefi et al.l |2008| ) , depending 
on the cooling processes of the intermediate muons, pions, and kaons. Recently, the use of fla- 
vor information has been especially proposed to extract some information on the particle physics 
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properties of the neutrinos and the properties of the source; see Ref . (jPakvasal 120081 ) for a review 



be used to extract information on the decay ( 


Beacom et al. 


2003; 


LiDari et al. 200' 


J; 


Maiumdari 


2007; 


Maltoni & Winter 2008; Bhattacharva et al. 


2009 


) and oscillation ( 


Farzan &: Smirnovll2002l; 


Beacom et al. 2004; 


Serpico & Kachelriefi 2005; 


Ser 


pico 


2006; 


3hattachariee &; Gupta 


2005: Winter 


2006; 


Maiumdar & Ghosal 


2007: Meloni h Ohlsson 


2007 


; Blum et al 


2007 


Rodeiohann 


2007: 


Xing 


2006; 


Pakvasa et al. 


2008: 


Hwang & Siveon 20081: Choubev et al. 2008: Esmaili Sz Farzai 


i 200J 


)) pa- 



rameters, in a way which might be synergistic to terre strial measurements . Of course, the flavor ra- 



tios m ay be also used for source identification, see, e.o.. lXing &; Zhoul (|2006l ); IChoubev k, Rodeiohann 



(|2009l ). Except from flavor identification, the differentiation between neutrinos and antineutri- 
nos could be useful for the discrimin ation between p7 and pp induced neutrino fluxes, or for the 
test of neutrino properties (see, e.g.. iMaltoni fc Winter! |2008| ) . A useful obs ervable may be the 
Glashow resonance process + e~ — > W~ — > anything at around 6.3 PeV ([Learned &; Pakvasa 



1995 



Anchordoqui et al 



2005 



Bhattachariee fc Guptall2005l ) to distinguish between neutrinos and 



antineutrinos in the detector. Within the photohadronic interactions, the tt + to tt~ ratio deter- 
mines the ratio between electron neutrinos and antineutrinos at the source. Therefore, a useful 
source model for these applications should include accurate enough flavor ratio predictions, includ- 
ing the possibility to include the cooling of the intermediate particles, as well as ir + to 7r~ ratio 
predictions. In addition, the computation of the neutrino fluxes should be efficient enough to allow 
for reasonable parameter studies or to be used as a fit model. 

Photohadronic interactions in astrophysical s ources are typ i cally e ither described by the refined 
Monte-Carlo simulation of the SOPHIA software iMiicke et al.l (120001 ) , which is partially based on 
Rachenl (119961 ). or are in very simplified approaches computed with the A-resonance approximation 



P + l 



n + 7r + 1/3 of all cases 
p + 7T° 2/3 of all cases 



(1) 



The SOPHIA software probably provides the best state-of-the-art implementation of the photo- 
meson production, including not only the A-resonance, but also higher resonances, multi-pion 
production, and direct (t-channel) production of pions. Kaon production is included as well by the 
simulation of the corresponding QCD processes (fragmentation of color strings). The treatment in 
Eq. ([T]), on the other hand, has been considered sufficient for many purposes, such as estimates 
for the neutrino fluxes. However, both of these approaches have disadvantages: The statistical 
Monte Carlo approach in SOPHIA is too slow for the efficient use in every step of modern time- 
dependent source simulations of AGNs and GRBs. The treatment in Eq. (TjQ), on the other hand, 
does not allow for predictions of the neutrino-antineutrino ratio and the shape of the secondary 
particle spectra, because higher resonances and other processes contribute significantly to these. 
One possibility to obtain a more accurate model is to use different interaction types with different 
(energy-dependent) cross sections and inelasticities (fractional energy loss of the initial nucleon). 
For example, one may define an interaction type for t he re s onances, and an interaction typ e for 



2TX 

multi-pion production (see, e.g., iReynoso &: Romero! ([2009J); iMannheim Sz Biermannl (11989) and 
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others). Typically, such approac hes do not distinguish 7r + f rom ir production. An interesting 



alternative has been proposed in iKelner &; Aharonianl (|2008l ). which approximates the SOPHIA 
treatment analytically. It provides a simple and efficient way to compute the electron, photon, 
and neutrino spectra, by integrating out the intermediate particles. Naturally, the cooling of the 
intermediate particles cannot be included in such an approach. 

Since we are also interested in the cooling of the intermediate particles (muons, pions, kaons), 
we follow a different direction. We propose a simplified model based on the very first physics 
principles, which means that the underlying interaction model can be easily adapted if new data 
are provided. In this approach, we include the intermediate particles explicitely. We illustrate the 
results for the neutrino spectra, but the extension to photons and electrons/positrons is straight- 
forward. 

More explicitly, the requirements for our interaction model are the following: 



The model should predict the 7r + , tt , and 7r° fluxes separately, which is needed for the 
prediction of photon, neutrino, and antineutrino fluxes, and their ratios. 

The model should be fast enough for time-dependent calculations and for systematic param- 
eter studies. This, of course, requires compromises. For example, compared to SOPHIA, our 
kinematics treatment will be much simpler. 

The particle physics properties should be transparent, easily adjustable and extendable. 

The model should be applicable to arbitrary proton and photon input spectra. 

The secondaries (pions, muons, kaons) should not be integrated out, because a) their syn- 
chrotron emission may contribute to observations, b) the muon (and pion/kaon) cooling affects 
the flavor ratios of neutrinos, and c) pion co oling may be in char g e of a second spectral break 
in the prompt GRB neutrino spectrum, see lWaxman fc Bahcalll (]1999i ). 



The cooling and escape timescales of the photohadronic interactions as well as the pro- 
ton/neutron re-injection rates should be provided, since these are needed for time-dependent 
and steady-state models. 

The kaons leading to high energy neutrinos should be roughly provided, si nce their different 
cooling properties may l e ad to changes in the ne utrino flavor ratios (see, e.o.. lKachelriefi &: Tomas 
20061 : iLipari et al.1 120071 : iKachelrieB et al.ll2008l ). Therefore, we incorporate a simplified kaon 
production treatment to allow for the test of the impact of such effects. This also serves as a 
test case for how to include new processes. 



For the underlying physics, we mostly follow similar principles as in SOPHIA. 

This work is organized as follows: In Sec. [2j we review the basic principles of the particle 
interactions, in particular, the necessary information (response function) to compute the meson 
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photo-production for arbitrary photon and proton input spectra. In Sec. O we summarize the 
meson photo-production as implemented in the SOPHIA software, but we simplify the kinematics 
treatment. In Sec. [H we define simplified models based on that. As a key component, we define 
appropriate interaction types such that we can factorize the two-dimensional response function. 
In Sec. we then summarize the weak decays and compare the different approaches with each 
other and with the output from the SOPHIA software. For the sake of completeness, we provide 
in App. [B] how the cooling and escape timescales, and how the neutron/proton re-injection can 
be computed in our simplified models. This study is supposed to be written for a broad target 
audience, including particle and astrophysicists. 



2. Basic principles of the particle interactions 



For the notation, we follow iLipari et al.l (|2007l ). where we first focus on weak decays, such as 
7r + — > fi + + Vp, and then extend the discussion to photohadronic interactions. Let us first consider 
a single decay process. The production rate Qb{Eb) (per energy interval) of daughter particles of 
species b and energy Ej, from the decay of the parent particle a can be written as 

/dn ^ 
dE a N a (E a ) T^ b (E a ) -^>( Ea ,E b ) . (2) 

Here N a (E a ) is the differential spectrum of parent particle^] (particles per energy interval), and 
r^\ b is the interaction or decay rate (probability per unit time and particle) for the process a — > b 
as a function of energy E a (which is assumed to be zero below the threshold). Since in pion pho- 
toproduction many interaction types contribute, we split the production probability in interaction 
types "IT". 

The function dn^ b / 'dE b {E a , E b ) describes the distribution (as a function of parent and daugh- 
ter energy) of daughter particles of type b per final state energy interval dE b . This function can 
be non-trivial. It contains the kinematics of the decay process, i.e., the energy distribution of the 
discussed daughter particle, other species, we are not interested in, are typically integrated out. If 
more than one daughter particle of the same species b is produced, or less than one (in average) 
because of other branchings, it must also give the number of daughter particles per event as a 
function of energy, which is often called "multiplicity" . 

Note that the decay can be calculated in different frames, such as the parent rest frame (PRF), 
in the center of mass frame of the parents (CMF), or in the shock rest frame (SRF), typically used 
to describe shock accelerated particles (such as from Fermi shock acceleration) in astrophysical 
environments. However, the cross sections, entering the interaction rate T, are often given in a 



In steady state models, this spectrum is typically obtained as the steady state spectrum including injection on 
the one hand side, and cooling/escape processes on the other side. 
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particular frame, which has to be properly included. In addition, note that Q and N are typically 
given per volume, but here this choice is arbitrary since it enters on both sides of Eq. ([2]). 

Sometimes it is useful to consider the interaction or decay chain a — > b — > c without being 
interested in b, for instance, for the decay c hain it —)•//—?• v. In this cas e, one can integrate out b. 



This approach has, for instance, been used in lKelner fc Aharonianl (120081 ) for obtaining the neutrino 



spectra. Note, however, that the parent particles a or b may lose energy before they interact or 
decay, such as by synchrotron radiation. We do not treat such energy losses in this study explicitely, 
but we provide a framework to include them, such as in a steady state or time-dependent model 
for each particle species. 

In meson photoproduction, a similar mechanism can be used. The interaction rate Eq. ([2]) can 
be interpreted in terms of the incident protons (or neutrons) because of the much higher energies in 
the SRF, i.e., the species a is identified with the proton or neutron, which we further on abbreviate 
with p or p' (for proton or neutron). In this case, the interaction rate depends on the interaction 
partner, the photon, as 



+1 

f f d cos 

r p7-> P '6(^p) = \ de \ r-^ 2 (1 - cos6» p7 )n 7 (e,cos6'p 7 )cr IT (e r .) . (3) 



-l 

Here n 7 (e, cos 6 P1 ) is the photon density as a function of photon energy e and the angle between 
the photon and proton momenta 9 P1 , cr IT (e r ) is the photon production cross section, and e r = 
E p e/m p {\ — cos p7 ) is the photon energy in the nucleon/parent rest frame (PRF) in the limit f3 p ~ 1. 
The interaction itself, and therefore E p and e, is described in the SRF. The daughter particles b 
are typically 7r + , 7r~, 7r°, or kaons. If intermediate resonances are produced, we integrate them out 
so that only pions (kaons) and protons or neutrons remain as the final states. 

In the following, we will assume isotropy n 7 (e, cos# P7 ) ~ n 7 (e) of the photon distribution. 
This limits this specific model to scenarios where we have seed photons produced in the shock 
rest frame (i.e., synchrotron emission). The other interesting case, where thermal photons are 
coming from outside the shock, may be easily implemented as well when the shock speed is high 
enough to assume delta-peaked angular distributions. Arbitrary angular distributions would require 
additional integrations. Assuming that the photon distribution is isotropic, the integral over cos # P7 
can be replaced by one over e r , we have 



oo 2E v e/m v 



p ^ pfh (E p ) = -^ / de^y- / de r e r a^(e r ). (4) 



P Hh™P 6th 
2E, } 



Here e t h — 150 MeV is the threshold below which the cross sections are zero. Note that e r corre- 
sponds to the available center of mass energy yfs of the interaction as 



s(e r ) = trip + 2m p e r 



(5) 
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In Eq. @, the integral over e r takes into account that the proton and photon may hit each other in 
different directions. In the most optimistic case, e r ~ 2E p e/m p (head-on hit), whereas in the most 
pessimistic case, e r ~ (photon and proton travel in same direction). It should be noted that the 
assumption of isotropy here is limiting the model to internally produced photon fields. This could 
be lifted for different angular distributions, but except for the case of uni-directional photons this 
would require one more integration. 

In general, the function dn^\ b /dEb(E a , Eb) in Eq. ([2]) is a non-trivial function of two variables 
E a and Eb- For photo-meson production in the SRF, the energy of the target photons is typically 
much smaller than the incident nucleon energy, and (3 a ~ 1. In this case0 

dn IT 

~^(E a , E b ) ~ 8(E b - X ^b E a ) ■ M b 1T . (6) 

The function X^a-^bi which depends on the kinematics of the process, describes which (mean) fraction 
of the parent energy is deposited in the daughter species. The function Ml describes how many 
daughter particles are produced at this energy in average. For our purposes, it will typically be 
a constant number which depends on interaction type and species b (if it changes at a certain 
threshold, one can define different interaction types below and above the threshold). For the 
relationship between xi^>b an d the inelasticity K (fractional energy loss of the initial nucleon) , see 
App. [Bj As we will discuss later, Eq. (0) is a over-simplified for more realistic kinematics, such as 
in direct or multi-pion production, since x a ^b i s > m g enerai > a more complicated function of Eb/E a 
and the initial state properties. In this case, the distribution is not sufficiently peaked around 
its mean, and the ^-distribution in Eq. © has to be replaced by a broader distribution function 
describing the distribution of the daughter energies. Instead of choosing a broader distribution 
function at the expense of efficiency, we will in these case define different interaction types with 
different values of X a ^-+b' simulating the broad energy distribution after the integration over the 
input spectra. 

As the next step, we insert Eqs. © and (HJ), valid for photoproduction of pions, into Eq. (|5J), 
in order to obtain: 

oo oo 

Qb(E b )= f^N p (E p ) [ de nj (e)R b (x,y) (7) 



E, 



E p 



with 



E h _E p e 
E p m p 



2 In fact, the ultra-relativistic argument justifies the introduction of a distribution function of only one final 
state variable dn%L> b /dE b {E a ,E h ) -> M b IT {E a )p{E b /E a ; X{ T (E a ), . . . , X{ T {E a )), where p is some parameterized 
probability distribution function of arbitrary shape and the parameters only depend on the initial state. The 
(^-approximation is the simplest approximation, which will only be useful if the probability distribution is sufficiently 
peaked around its mean. 
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and the "response function" 

R b ( x ,y)^^2R^( x , y )^Y,Y2 I de r e r a 1T (e r )Ml T (e r )5(x- X 1T (er)) • (9) 

IT IT ^ J 

Here < x = Eb/E p < 1 is the secondary energy as a fraction of the incident proton or neutron 
energy, eth < 2y = 2E p e/m p < 10 4 GeV corresponds to the maximal available center of mass 
energ)|f|, and the x IT is the fraction of proton energy deposited in the secondary. Note that x IT ( e r) 
and Mb(e r ) are typically functions of e r , if they only depend on the center of mass energy of the 
interaction. We will define our interaction types such that Mj is independent of e r . 

If the response function in Eq. @ is known, the secondary spectra can b e calculated from 



Eq. ([71) for arbitrary injection and photon spectra. A similar approach was used in lKelner &i Aharonian 



(|2008l ) for the neutrinos and photons directly. However, we do not integrate out the intermediate 
particles, which in fact leads to a rather complicated function R(x,y), summed over various inter- 
action types. We show in Fig. Q] the cross section as a function of e r for these interaction types 
separately, where the baryon resonances have been summed up. Naively, one would just choose the 
dominating contributions in the respective energy ranges in order to obtain a good approximation 
for <7 IT (e r ). However, the different contributions have different characteristics, such as different 7r + 
to tt~ ratios in the final states, and therefore different neutrino-antineutrino ratios. In addition, 
the function 5(x — x IT ( e r)) in Eq. ([9]) maps the same region in e r in different regions of x = Ef,/E p , 
depending on the interaction type. This means that while a particular cross section may dominate 
for certain energies e r , for instance, the pion energies where the interaction products are found 
can be very different, leading to distinctive features. Therefore, for our purposes (such as flavor 
ratio and neutrino-antineutrino predictions), it is not a sufficient approximation to just choose the 
dominating cross section. 

From the particle physics point of view, Rb(x,y) = Rb(Eb, E p ,e) is very similar to a detector 
response function in a fixed target experiment. It describes the reconstructed particle energy 
spectrum as a function of the properties of the incident proton beam. As the major difference, the 
"detector material" is kept as a variable function of energy, leading to the second integral over the 
photon density. 



3. Review of the photohadronic interaction model 



The description of photohadronic interactions is based on iRachenl (|l996l ) ; iMucke et al.l (|2000l ) , 
which means that the fundamental physics is similar to SOPHIA. However, our kinematics will be 
somewhat simplified. The purpose of this section is to show the key features of the interaction 



3 This range is given by the range the cross sections are known. For higher energies, our model can only be used 
as an estimate. 
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1000 rr 




e r [GeV] 



Figure 1 The total p7 photo-meson cross section as a function of the photon's energy in the proton rest frame 
e r analog to iMiicke et al, (1/ibarn = 10 30 cm 2 ; data, shown as dots, from [Particle Data Groupl (|2008h )V 

The contributions of baryon resonances (red, dashed), the direct channel (green, dotted) and multi-pion production 
(brown) are shown separately. 



model. In addition, it is the first step towards an analytical description for the response function 
in Eq. ([9]), which should be as simple as possible for our purposes, and as accurate as necessary. 
Note that we do not distinguish between protons and neutrons for the cross sections (there are 
some small differences in the resonances and the multi-pion production). 



3.1. Summary of processes 

In summary, we consider the following processes: 

A-resonance region The dominant resonance process is the A(1232)-resonance (at e r = 340 MeV; 
c/., Eq. © for the e r -s-conversion): 

P + 7 — > P + 7T • (10) 

This process produces neutral (for p' = p) or charged (for p' 7^ p) pions. For instance, for 
protons in the initial state, see Eq. (UJ. 
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Higher resonances The most important higher resonance contribution is the decay chain 



7 + p^i V A / + vr, A'^p' + vr' (11) 
via A- and ^-resonances. Other contributions, we consider, come from the decay chain 

1+P -> P+P > P^TT + TT ■ (12) 



Direct production The same interactions as in Eq. (|10p or Eq. (jlip (with the same initial and 
hnal states) can also take place in the t-channel, meaning that the initial 7 and nucleon 
exchange a pion instead of creating a (virtual) baryon resonance in the s-channel, which 
again decays. This mechanism is also called "direct production", because the properties of 
the pion are already determined at the nucleon vertex. For instance, the process p+7 — > p+7r° 
only takes place in the s-channel, whereas in the process p + 7 — > n + ir + the t-channel is 
possible as well, because the photon can only couple to charged pions. The s-channel reactions 
typically have a pronounced peak in the e r -distribution, whereas they are almost structureless 
in the momentum transfer distribution. On the other hand, t-channel reactions do not have 
the strong peak in the center of mass energy, but have a strong correlation between the initial 
and final state momentum distributions, implying that the pions are found at lower energies. 
On a logarithmic scale, however, the momentum distribution of the pions is broad. 

Multi-pion production At high energies the dominant channel is statistical multi-pi on produc- 
tion leading to two or more pions. The process is in SOPHIA iMiicke et al.1 (|2000l ) described 
by the QCD-fragmentation of color strings. 



The effect of kaon decays is usually small. However, kaon decays may have interesting conse- 
quences for the neutrino flavor ratios at very high energies, in particular, if strong magnetic 
fields are present (the kaons are less sensitive to sync hrotron losses because of the larger rest 
mass) (jKachelriefj &: Tomasll2006l : iKachelriefi et al.ll2008l ). Therefore, we consider the leading mode: 
K + production (for protons in the initial state) with the decay channel leading to highest energy 
neutrinos@ Note that at even higher energies, other processes, such as charmed meson production, 
may contribute as well. 

We show i n Fig. [J the total jry photo- meson cross section as a function of e r , analogously 
Miicke et all (|200cl ) (1/nbarn = 1CT 30 cm 2 ; data from IParticle Data Group! (|2008l ))). The con- 



to 

tributions of baryon resonances (red, dashed), the direct channel (green, dotted) and multi-pion 
production (brown) are shown separately. 

In order to fully describe Eq. Q for each interaction type, we need the kinematics, entering 
X IT , the multiplicities M^ T , and the cross section <t it . 



4 The contributions from K and K° are about a factor of two suppressed iLipari et al.l ( 20071 ') . and K + has a 
leading two body decay mode into neutrinos. 
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3.2. Kinematics and secondary multiplicities 



The kinematics for the resonances and direct production can b e effect ively described in the two- 
body picture. We follow the calculation of iBerezinskv fc Gazizovl (|1993l ) for the reaction p + 7 — > 



p' + 7r with E p 3> Ey (/3p ~ 1) in the SRF. The Lorentz factor between CMF and SRF is: 



E r> + 



E 



7c 



(13) 



with sfs the total CMF energy from Eq. ((5|) . The energy of the pion in the SRF can then be written 
as: 



7cm£r< 



-1 1 ocm ncn 
1 + p n COS tin 



En 



EE 



(l + ^cosO 



(14) 



with E™, and #£. m the pion energy, momentum and angle of emission in the CMF. Note that 
backward scattering of the pion means that cos^ m = —1. From Eq. (|14p we can read off the 
fraction of energy of the initial p going into the pion E^/Ep (which is equal to the inelasticity of 
p f or p' = p). Since p + 7 — > p ' + tt is a two body reaction, the energies of p' and pion are given 



bv lParticle Data Groupl ()2008l ) 



El 



nip + m\ 



2ji 



with s(e r ) given by Eq. ([5]). Using Eq. (Tl5l) in Eq. (fl4|) . we can calculate x i n Eq. 

s(e r ) 



as 



m,p + m. 2 



2s(e r ) 



^(l + ^cos^) 



(15) 



(16) 



Indeed, x( e r) is a function of e r if cosf?,,- is a constant. As described in iRachenl (| 19961 ). the average 
(cos^tt) ~ for the resonances (such as for isotropic emission in the CMF), whereas the direct 
production is backward peaked to a first approximation for sufficiently high energies (cos 6*^) — > 
— 1. This approximation is only very crude. Therefore, we calculate the mean (cos^) for direct 
production by averaging the probability distribution of the Mandelstam variable t, as we discuss in 
App.|A] Note that even for the resonances these scattering angle averages are only approximations; 
a more refined kinematics treatment, such as in SOPHIA, will lead to a smearing around these mean 
values. In addition, note that we also use Eq. (I16p for the kaon production (with the rough estimate 
(cos (9^) ~ 0), where the replacements m p —> m\ and —> ra u have to be m ade. Similarly, one 
has for the first and second pion for the interaction in Eq. (IllD (|Rachenlll996l ) 



Xa(e r ) 

Xb(tr) 



s(e r ) 



m A + m 



2s(e r ) 
s(e r ) -ml + m\ m\ 



^(l + ^cos^) 



m 2 + m 2 



2s(e r ) 

where A = A(1232), and for the interaction Eq. ([12 



2m\ 



-(I + ^cos^a) 



1 s(e r ) — rnl + rn 2 



2s(e r ) 



p{ i + ^ m cose. 



K 

T 



(17) 
(18) 

(19) 
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Table 1 Summary of the average multiplicities for different types of considered resonances and direct production 
channels. The multiplicities of the pions M n add up to the number of pions produced (one or two), the multiplicities 
for the nucleons M' v in the final state p to one. In the last columns, the cosine of the approximate scattering angle 
in the center of mass frame is given, together with the equation for kinematics. The interaction types IT are labeled 
by the type of resonance (A or iV) or "Direct" for the direct production (f-channel). 

with m p ~ 775 MeV. Note that in order to evaluate the (^-function in Eq. ([9]) (if one integrates over 
e r ), one also needs the derivative x'( e r), which can easily be computed from the functions above. 

The average multiplicities for different types of resonances and direct production channels can 
be found in Table [TJ The multiplicities of the pions M n add up to the number of pions produced 
(one or two), the multiplicities for the nucleons Mp in the final state p' to one. They are needed in 
Eq. ([6]) and Eq. (IBlj) (cooling and escape timescale in App. [B]). In the last columns, the cosine of 
the approximate scattering angle in the center of mass frame is given, together with the equation 
for kinematics. The interaction types IT are labeled by the type of resonance (A or N) or T for the 
direct production (t-channel). Note that different resonances will contribute with a similar pattern, 
such as through the interaction types Ai(1232), Ai(1700), Ai(1905), Ai(1950), etc.. In addition 
note that the branching ratios into certain resonances and from a resonance into the described 
channels are absorbed in the cross sections. Therefore, the total yield is always one (or two) pions 
per interaction in Table [TJ 



For the multi-pi on channel, we use two d ifferent approaches. A very simple treatme nt can be 
performed similar to lAtoyan &: Dermerl (120031 ). with some elements of lMiicke et al.l (J200CJ). Most of 



-13- 



the energy lost by the proton (~ 0.6 E p , or K ~ 0.6) is split equally among three pions. The three 
types 7r + , 7T~, and ir° are therefore approximately produced in equal numbers. Our multi-pion 
channel, parameterized by the multi-pion cross section in Fig. [TJ however, is actually a combination 
of different processes and residual cross sections. For instance, diffractive scattering is a small 
contribut ion. Therefor e , in or der to reproduce the pion multiplicities times cross sections in Figs. 9 
and 10 of lMiicke et al.1 ([2000) more accurately, we choose M w o = 1 (M n o = 1), M w + = 1.2 (M w + = 
0.85), and M n - = 0.85 (M n - = 1.2) for initial protons (neutrons). Close to the threshold, the 
decreasing phase space for pion production requires a modification of the threshold e r > 1 GeV for 
7r~ (vr + ) from initial protons (neutrons). This corresponds to a vanishing multiplicity below the 
threshold, i.e., we assume that below e r > 1 GeV only two pions are produced. We assume for the 
fraction of the proton energy going into one pion produced in the multi-pion channel x Multi -' r ~ 0.2. 
In addition, we choose M p = 0.69 and M n = 0. 31 for initial proton s (or M p = 0.31 and M n = 0.69 for 



initial neutrons) in accordance with Fig. 11 of lMiicke et al.1 (|2000i ) for high energies. As alternative, 
we use a more accurate but computationally more expensive approach, which is directly based on 
the kinematics of the fragmentation code used by SOPHIA (c/., Sec. 14.4.2]) . 



3.3. Cross sections 



We parameterize the cross sections of photohadronic interactions following iMiicke et al.1 (J2000J). 
We split the processes into three parts: resonant, direct, and multi-pion production. The different 
contributions are shown in Fig. [TJ 

The low energy part of the total cross section is dominated by excitations and decays of baryon 
resonances N and A. The cross sections for these processes can be described by the Breit-Wigner 
formula 



cr Bw( £ r) 



m 



2\2 



4tt(2J + l)B y B out sT 2 
(s-M 2 ) 2 + sT 2 



B 



IT 

° Ut e2 (s 



( M IT)2)2 + ( r lT)2 g 



(20) 



with s(e r ) given by Eq. ([5]). Here J, M, and T being the spin, the nominal mass, and the width 
of the resonance, respectively. We consi der the energy-depe ndent branching ratios B OVLt and the 
resonances shown in the Appendix B of IMiicke et al.1 (|2000l ) ; the total branching ratios are also 
listed in Table [2j Note that the energy-dependent branching ratios respect the different energy 
thresholds for different channels (e.g., interaction types 1 to 4). For simplicity, we neglect the 
slight cross section differences between nj- and /^-interactions and use the values for protons, 
which implies that we do not take into account the N(1675) resonanceJ§ In Table [2] we list all 
considered resonances and their contributions to the different interaction types. To account for the 



For a more accurate treatment for neutrons, also include the N(1675) resonance, whereas N(1650) and N(1680) 
do not apply. In addition, in Eq. ((25)) . replace 80.3 ->■ 60.2 and in Eq. J26), replace 29.3 — > 26.4. 
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Properties 


Z 


Total B ont 


Resonance 


M [GeV] 


T [GeV] <7 


[/ibarn] 


e th [GeV] w 


[GeV] 


IT 1 


IT 2 


IT 3 
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1.231 


0.11 
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0.152 


0.17 
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1.440 


0.35 


1.389 


0.152 


0.38 


67% 


33% 




JVi T (1520) 


1.515 


0.11 


25.567 
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0.38 


52% 


42% 


6% 


AT IT (1535) 


1.525 


0.10 


6.948 


0.152 


0.38 


45% 






ATi T (1650) 


1.675 


0.16 


2.779 


0.152 


0.38 


75% 


14% 




ATit(1680) 


1.680 


0.125 


17.508 


0.152 


0.38 


64% 


22% 


14% 


Ait (1700) 


1.690 


0.29 


11.116 


0.152 


0.38 


14% 


55% 


31% 


Ait(1905) 


1.895 


0.35 


1.667 


0.152 


0.38 


14% 


13% 


73% 


A IT (1950) 


1.950 


0.30 


11.116 


0.152 


0.38 


37% 


39% 


24% 



Table 2 Summary of the considered resonances. Properties of the resonances and parameters for the quenching 
function Eq. ((21} taken from Miicke et al. (2000;). The four rightmost columns refer to possible interaction types IT 
in Table [T] and show the total branching ratios. In fact, at the end, each combination of a specific resonance and a 
particular decay chain corresponds to a separate interaction type, such as Ai(1232). 



phase-space reduction near threshold the function Z lT is introduced and multiplied on Eq. (i20j) : 



IT TTn 



if e r < e IT 



th ' 



1 



ifeJi <e r <w"+e£, 
if e r >u; IT + e[ T . 



(21) 



with ej£ and ui IT taking the values shown in Table [2j 

For the direct and the multi-pion cross section we adopt the formulae from SOPHIA. The 
direct part is given by 

(e r -0.29) 2 \ / (e r -0.37) 2 ' 



ff Al (er) = 92.7Pl(e r , 0.152, 0.25, 2) +40exp 
a T2 (e r ) = 37.7Pl(e r ,0.4,0.6,2) 
with 

Pl( e r> £th, tmax, «) : 



0.002 



15 exp 



0.002 



( 6r ~ 



A-a 



if e r < e t h 
else 



(22) 
(23) 

(24) 



where A = a.e max /€^- The cross sections are given in //barns, and the interaction types are taken 
from Table [H i.e., Ti is the single pion direct production, and T2 is the two pion direct production. 
The multi-pion cross section can be parameterized by the sum of the following two formulas (cross 
sections in //barns): 



_Multi — 7T1 \ 



a 



Multl-7T2 , 



e r ) 



80.3Q / (e r ,0.5,0.1)s~ 
e r - 0.85 



0.31 



1 — exp 



0.69 



(29.3 s" - 34 + 59.3 s - 095 ) if e r > 0.85. 



(25) 
(26) 
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We have compared our resonant and direct production multiplicities times cross sections and 
proton to neutron ratios with SOPHIA (c/., Figs. 9 to 11 in lMiicke et al.l (|2000l )). and they are in 
excellent agreement. 



4. Simplified models 

Here we discuss simplified models based on Sec. [3l where we demand the features given in the 
introduction. For the resonances, we propose two methods, one is supposed to be more accurate, 
the other one faster. 



4.1. Factorized response function 

First of all, it turns out to be useful to simplify the response function in Eq. Q. In our 
simplified approaches, we follow the following general idea: In Eq. ([7]) and Eq. ([9]), in principle, (at 
least) two integrals have to be evaluated. Let us now split up the interactions in interaction types 
such that x IT ( e r) = X IT an d M^ T (e r ) = M^ T are approximately constants independent of e r , but 
dependent on the interaction type. The response function in Eq. (J9]) then factorizes in 

2y 

R 1T (x, y) = 5(x - x 1T ) Mf f lT (y) with f lT (y) = ± J de r e r a 1T (e r ) . (27) 

The part 5(x — x IT ) describes at what energy the secondary particle is found, whereas the part 
M^ T / IT (y) describes the production rate of the specific species b as a function of y. If only 
the number of produced secondary particles is important, it is often useful to show the effective 
F b (y)^EiTMl T f 1T (y). 

As the next steps, we evaluate Eq. (J7J) with Eq. (p7|) by integrating over dx/x = —dE p /E p and 
by re- writing the e-integral as y = E p e/m p -mtegral: 

oo 

QV = N p ^ J dyn, (^f^) M fe IT / IT ( y ) . (28) 

W 2 

This (single) integral is relatively simple and fast to compute if the simplified response function 
/ IT (y) together with x IT is given. Therefore, the original double integral simplifies in this single 
integral, summed over a number of appropriate interaction types. In the following, we therefore 
define M^ T , / IT (y) and x IT for suitable interaction types. 
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4.2. Resonances 

Here we describe two alternatives to include the resonances. Model A is more accurate but 
slower to evaluate, because it includes more interaction types, whereas model B is faster for com- 
putations, since there are only two interaction types defined for the resonances. 



4-2.1. Simplified model A (Sim- A) 

In this simplified resonance treatment we make use of the fact that the resonances have cross 
sections peaked in e r . In principle, these can be approximated from Eq. (|20f) and Eq. (|2ip by a 
5-function 

oo 

a IT (e r ) ~ BZ f it 5(e r - e* T '°) with f IT = J <jg w Qf(e r , e\l, w 1T )de r , (29) 

o 

where f IT and er T ' correspond to surface area (in ^barn-GeV) and position (in GeV) of the 
resonance, and are process-dependent constants. Therefore, using Eq. (|29p . the e r -integral in Eq. ([9]) 
can be easily performed, and we obtain again the simplified Eq. (|27p with 

/ IT (2/) = ^4 T '°^tf IT e(2 y -4 T '°) and X IT = ^it.O) (30) 

The relevant parameters for all interaction types can be read off from Table El Note that M* T are 
the total (i.e., energy-independent) branching ratios from Table [2l In addition, note that in some 
cases, the resonance peak may be below threshold for some of the interaction types. However, for 
reasonably broad energy distributions to be folded with and the total branching ratios used, this 
should be a good approximation. The pion spectra are finally obtained from Eq. (|28j) ■ summing 
over all interactions listed in Table El 



4-2.2. Simplified model B (Sim-B) 

In this model, we take into account the width of the resonances and approximate them 
by constant cross section s within certain energy ranges. Compared to approaches such as in 



Revnoso fc Romero! (|2009l ). we distinguish between ir + and it production (i.e., do not add these 



fluxes) and include the effects from the higher resonances. 

From Fig. [TJ one can read off that there are two classes of resonances: The first peak in Fig. [I] is 
dominated by the A(1232)-resonance (lower resonance - LR), whereas the higher resonances (HR) 
contribute at larger energies. In addition, the kinematics (c/., x m Table EJ) and the multiplicities 
(c/., Table [1]) are very different. For example, protons and neutrons interactions via the A(1232)- 
resonance happen through only one interaction type, and produce either only tt + or only tt~ (and 
7r ). For the higher resonances, the pions are produced in all pion charges. 
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Table 3 Parameters for the resonant pion production in our simplified treatment A. The units of e' T '° are GeV, 
the units of f IT are ^tbarn GeV. 

For the resonances, we define two interaction types LR and HR, as shown in Fig. [2] (boxes). 
The interaction type LR corresponds to Ai(1232), whereas the interaction type HR contains the 
higher resonances. The properties of these interaction types are summarized in Tabled! 

The surface area covered by these cross sections, corresponding to the T in simplified model A, 
is chosen for LR and HR such that the pion spectra match to these of the previous section for 
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Figure 2 Cross section for the resonances as a function of e r (thick curve). Green (light gray) resonances are 
A-resonances, blue (dark gray resonances) are iV-resonances. The total resonance cross section is shown as thick 
curve. Our simplified model B is depicted by the red (gray) boxes, where the interaction types are labeled. The 
dashed curve refers to our simplified cross section for kaon production "KP" . 

typical power law spectra (with spectral indices of about two). Of course, there can never be 
an exact match, because the contributions from the individual resonances depend on the spectral 
index. However, as we will demonstrate, this estimate is good enough for our purposes. 

The averaged numbers for the multiplicities and inelasticities for the higher resonances interac- 
tion type are estimated from the interaction rate Eq. (j4]) by assuming that all resonances contribute 
simultaneously and weighting by the interaction type-dependent part er T '° f IT B^ t (c/., Eq. (1B3P 
using Eq. (j30l) : c/., App.[B]). By the same procedure we obtain the average x~ values from Eq. ([9]). 
Note that the tt~ (for proton interactions) are, in average, reconstructed at lower energies, because 
these are mostly produced by interaction type 2, which has smaller x-values (c/., Table [3|). In 
addition, note that the total number of charged pions is, in fact, close to one per interaction here, 
and the total number of pions close to 1.5, since in interaction types 2 and 3 more than one pion 
is produced. 
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The pion spectra are then computed from Eq. ()28p with 
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(31) 



(32) 



as it can be shown from Eq. (|27p . Here M^ T and x IT (needed in Eq. ()28p ) are chosen according to 
interaction type, initial proton or neutron, and final pion, as given in Table [H 



4-2.3. Comparison of resonance response functions 

In Fig. O we compare the response function F 7T (y) = Xjt f lT {u) between simplified model 
Sim- A (thin solid curves) and Sim-B (thick curves), summed over the resonances. This function is 
proportional to the number of produced pions of a certain species as a function of y, whereas the 
rr-dependent part in Eq. (|27p describes the energy at which the pions are found. Obviously, the 
function is much smoother for Sim-B than for Sim-A. Because of only a few contributing interaction 
types and the smoothness of the function, the evaluation will be much faster. Once the photon 
spectrum is folded in, the contributions of both response functions will be very similar. 

Note that model Sim-B includes the part of the A(1232)-resonance below the peak. This can 
be seen by comparing with the dashed curve, which represents model Sim-A but the interaction 
type Ai(1232) replaced by the full Breit-Wigner-form. 
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Table 4 Parameters for the resonant pion production in our simplified treatment B. 
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Figure 3 Response function F n (y) = ^2 iT f IT (y) summed over the resonances only. Here simplified model 
Sim- A (thin solid curves) is compared with simplified model Sim-B (thick curves). The thin dashed curves correspond 
to Sim-A with the full Breit-Wigner form of the interaction type Ai(1232) only. 
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Table 5 Parameters for the direct pion production in our simplified treatment. The inelasticity for T2b is included 
in T 2a . 



4.3. Direct production 

For direct production, we also follow the approach in Sec. 14.11 However, this interaction type 
is tricky, since the kinematics function x IT ( e r) is strongly dependent on the interaction energy (see 
App.|A]for details). This implies that the ^-distribution in Eq. ([6]) is not a good approximation and 
needs to be replaced by a broader distribution function. In addition, the distribution of scattering 
angles will lead to smearing effects. In this case, it turns out to be a good approximation to split 
the direct production in different interaction types with different characteristic values of % IT as a 
function of e r , which simulates such a broad distribution after the integration of the input spectra. 
We define three interaction types Til, ^Ym, and Tin for direct one pion production and four for 
two pion production, namely T2 a L, T^aM, and T2 a H for the first pion, and for the second pion. 
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The interaction types are shown in Table and the names correspond to Table Q] split into low 
(L), medium (M), and high (H) energy parts in e r , respectively, limited by ej^ in and e|^ ax . After 
this splitting, the additional effect of the scattering angle smearing turns out to be small, as we 
will demonstrate later. We obtain the function / IT for these interaction types from Eq. (|27p as 



IT 





1 
1 

W 1 



(/ IT (4 T ax) 



with / IT (2y) re-parameterized using x = log 10 (y/GeV): 



2y <e 
e IT < 

mm — 



' IT (4 T iJ) 2y>6 



IT 

oin 

2y 

IT 



< e 



IT 



(33) 



/ T2 (2y) 





2y < 0.17 GeV 
35.9533 + 84.0859x + 110.765x 2 + 102.728x 3 + 40.4699a; 4 

0.17 GeV <2y < 0.96 GeV 
30.2004 + 40.5478x + 2.03074a; 2 - 0.387884x 3 + 0.025044x 4 

2y > 0.96 GeV 

2y < 0.4 GeV 

-3.4083 + + 40.7160 In (2y) 2y > 0.4 GeV 



(34) 



(35) 



Note that e^ in , e^, the multiplicities and x lT can be read off from Table [5] for the different inter- 
action types. The integral values I 1T are the same for interaction types Til, Tim, Tih (Eq. (|34"|) ). 
and for interaction types T2 a L> T2 a M 5 T2 a H 5 and T 2 b (Eq. ([35]) ). In Eq. ([28]) . all the interaction 
types in Table [5] have to be summed over. 



4.4. Multi-pion production 

Here we show two different approaches to the multi-pion channel. The first, simplified approach 
will later be shown as model Sim-C, the second, more refined approach will be used in all other 
models. 



4-4- 1- Simplified kinematics 

We start with the simplest example for multi-pion production, for which we follow Sec. 13.21 We 
assume X M ^-* ~ 0.2, i.e., the response function factorizes as R u ^^ = 5(x-0.2) M™ nm ~* / Multi-7r (y). 
We obtain from Eq. fl27|) and Eq. ([28]) 



oo 

QMuIti-^) = ^(5^) J dy M MuM-. /Multi-^) ^ (36) 

eth/2 
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with 

22/ 

f UxM -*(y) = ^Jde r tr <T MuM -*(e r ) . (37) 

As described in Sec. H we assume that M w o = 1 (M^o = 1), M 7r + = 1.2 (M T+ = 0.85), and 
M n - = 0.85 (M„.- = 1.2) for initial protons (neutrons). The function / Multi_7r (y) can be obtained 
using the sum of Eq. ([250 and Eq. ([260 . We numerically integrate it and re-parameterize it with 
x = log 10 (y/GeV) by 

2y < 0.5 GcV 
0.5 GcV < 2y < 14 GcV 
2y > 14 GeV 

(38) 



f{x) v - I 73.9037 + 187.526a; - 161.587a; 2 - 206.268a; 3 + 
^ibarn ~~ | +354.02a; 4 - 129.759a; 5 

[ 131.839 - 25.3296a; + 10.612a; 2 - 0.858307a; 3 + 0.0493614a; 4 

for 7r~ (initial proton) or ir + (initial neutron). Note that the different function for ir~ comes from 
the different threshold below which we have set the cross section to zero, as discussed in Sec. 13.21 
For y 3> 10 4 GeV, this function still increases as extrapolation of Fig. Q] (the cross section is not 
measured in that energy range). 





87.5538 + 120.894a; - 98.4187a: 2 - 59.6965a; 3 + 67.2251a; 4 
131.839 - 25.3296a; + 10.612a; 2 - 0.858307a; 3 + 0.0493614a; 4 



/xbarn 

for 7T°, and tt + (initial proton) or 7r~ (initial neutron), and 



2y < 1 GeV 

1 GcV < 2y < 10 GeV 
2y > 10 GeV 



(39) 



4-4-%- Kinematics simulated from SOPHIA 

Similar to the direct production, the (5-distribution in Eq. (jSJ) is not a good approximation 
and needs to be replaced by a broader distribution function. We again solve this by splitting 
the multi-pion production in different interaction types with different characteristic values of x IT > 
which simulates such a broad distribution after the integration of the input spectra. Compared to 
the direct production, the cross section can be parameterized relatively easily by step functions. 
However, the pion multiplicities in fact increase with energy. This means that the higher the 
interaction energy, the more pions are produced, which (in average) are found at lower energies. In 
addition, the final energy distribution functions are broad. 

We solve this strongly scale-dependent behavior by dividing the e r -range in seven interaction 
types Mj, each with a particular average cross section and average pion multiplicities, which we 
directly take from SOPHIA. In addition, we split the pions in each of these samples in a lower energy 
(L) and higher energy (H) part, which are reconstructed at different values of x IT to simulate the 
broadth of the distributions within the same e r . Typically, the L-sample corresponds to the peak 
of the distribution (at least for high energies), whereas the H-sample simulates the tail of the 
distribution. The splitting of the multiplicities can be performed automatically once a splitting 
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Table 6 Parameters for the multi-pion production in our SOPHIA-based treatment. The inelasticities for the 
H-sample are included in the L-sample. 



point is denned. We choose the x _vaiues °f the interaction types to simulate the peaks of the 
distribution and to reproduce the total energy going into pions in SOPHIA. This treatment is, of 
course, not extremely accurate, but it allows to use the fast approach in Eq. (|27p while obtaining 
accuracy at high energies. 



The function / M can be easily calculated from Eq. 
be constant within each IT: 



fM, 







(2j/) 2 



\iy) 2 - (^n) 2 ) 

V^-max/ \ mm/ / 



271) since the cross section is assumed to 



(40) 



< 2y < e max 

2y > i 



The interaction types are listed in Table El where the parameters for this equation can be found, 
as well as the x an d K values and multiplicities for the individual interaction types. 



4.5. Kaon production 

Here we include K + production into our simplified model. If needed, other (sub-leading) 
channels can be added to our framework as well (such as K~ production if the neutrino-antineutrino 
ratio in the higher energy regime is studied). This example should also serve as illustration how to 
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add processes to our model. Note that this implementation, however, is more primitive than the 
pion production above. 

The production of K + in photohadronic interactions cannot be assigned to a single resonance, 
but comes from a number of resonances mainly decaying into K + and A or E° (with relatively sim- 
ilar masses, which means that we can neglect the mass difference). In addition, at high interaction 
energies, multi-fragmentation, similar to multi-pion produc tion, contributes significant ly. The pro- 
duction cross section up to e r ~ 2 GeV has been measured (jSaphir Collaborationlll998l). For higher 



energies, where no data are available, one may use extrapolations f rom models (|Lee et al.l 12001 



to 



Asano fc Nagatakl 2006 ) . Approximating the data in Lee et al. ( 2001 ) and extrapolating according 



Lee et al.l (|200ll ). we define an interaction type "KP" and model the total K + production cross 



section as (c/., dashed curve in Fig. [2]) 



KP 



e r <1.0GeV 

2.0 //barn 1.0 GeV < e r < 1.2 GeV 

3.7//barn 1.2 GeV < e r < 1.65 GeV 

2.7 /ibarn e r > 1.65 GeV 

Accordingly, we obtain the response function from Eq. (|27p with M K + = 1: 



(41) 



pKP 





2.0 /ibarn 
3.7 /ibarn 
2.7 /ibarn 



1.0 GeV 2 

0.3 GeV 2 ' 

(2j/) 2 
0.16 GeV 2 ' 

(2y) a , 



2y < 1.0 GeV 

1.0 GeV < 2y < 1.2 GeV 

1.2 GeV < 2y < 1.65 GeV 

2y > 1.65 GeV 



(42) 



From Sec. 13.21 we have Xk+ — 0-35 (computed for e r ~ 1.4 GeV at the peak). The effects on primary 
cooling and secondary re-injection in App. [B]are negligible. Note that the absolute normalization 
of the kaons at high energies crucially depends on the extrapolation of the cross sect ions to high 



energi es. Our kaon production is about a factor of two smaller at higher energies than in lLipari et al 
(|2007l ). since we assume the cross section to be constant at high energies. 



4.6. Comparison of individual contributions 

In Fig. HI we show the individual contributions to F w (y) = M^ T / IT (y) for pion production 
as a function of y, which is proportional to the maximal available center-of-mass energy. This 
quantity describes the total number of pions of a particular species produced as a function of y, 
including multiplicity and cross section. However, it does not include the energy where the pions 
are found. For initial protons, tt + are most abundantly produced because of the direct production 
contribution, followed by ir° and then tt~ . For 7r + production close to the threshold, the direct 
production and resonances are most important, whereas for 7r° production, the direct production 
hardly has any impact, and the resonances dominate. For ir~ production, at the threshold all 
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Figure 4 The function F n (y) = J2it M ^ f T (v) for 

pion production (resonance treatment from simplified model B), 

with the individual contributions marked. 



processes contribute almost equally, only the lower resonance does not take part. For larger y, in 
all cases multi-pion production dominates. 

As far as the relative contribution is concerned, thanks to the direct production, the ir + are 
always produced most abundantly. As one can read off from Fig. [H this statement is independent 
of the photon or proton spectra used, because the response function for tt + is always larger than 
the one for 7r°, in a large part of the energy range even about 50% larger. This means that for 
arbitrary proton and photon spectra, thanks to the direct production, the tt + :tt ratio lies between 
about 1:1 and 3:2, whereas the A(1232)-approximation in Eq. ([1]) predicts 1:2. In fact, one can 
show that the minimum of the charged to neutral pion ratio with respect to y is 

mm — — — ~ 1.2 (43) 

v FAV) 



for arbitrary input spectra; see Miicke et al. (200C) for a specific photon field |§ From this equa- 



tion, any neutrino flux computed with the A(1232)-approximation and normalized to the observed 
photon spectrum, if mainly coming from 7r° decays, is underestimated by a factor of at least 
1.2/0.5 = 2.4, i.e., the neutrino flux should be about a factor of 2.4 larger. Even the often used ap- 
proximation that 50% of all photohadronic interactions result in charged pions underestimates the 
neutrino flux by at least 20%. These numbers are to be interpreted as lower limits for the neutrino 
flux underestimation, the exact values depend on the input spectra and may be even higher. 

Of course, the overall impact of the individual contributions depends on the proton and photon 
spectra the response function is to be folded with. In order to check the impact of the individual 



6 As mentioned above, F 1T (y) does not include the reconstructed energy of the pions. Since, however, the pions of 
the different species are found in similar energy ranges, this statement also roughly applies to the energy deposited 
into the different species. 
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contributions on typical AGN, GRB or high temperature black body (BB) spectra, we define three 
benchmark spectra in App. [Cj Note that all of the spectra shown in this work are given in the 
SRF. In Fig. we show the ir + , ir~, and tt° spectra for these benchmarks. The upper panels 
are for the GRB benchmark, the middle ones for the AGN benchmark and the lower ones for the 
BB benchmark. For the GRB benchmark, direct production dominates at low energies for the 7r + 
production, whereas the 7r° production at low energies is determined by the resonances. The ir~ 
production is dominated at all energies by the multi-pion production, such as the other spectra at 
high energies. Therefore, all different processes are important, but they contribute entirely different 
to the different pion polarities. One can also read off from this figure, that the characteristic shape 
of the GRB pion or neutrino spectra often shown in the literature (resonance curves), which is flat in 
the middle energy range, is tilted upwards due to multi-pion production. For the AGN benchmark 
(middle panels), the tt° production is governed by the resonances in the whole energy range, the 
7r~ production by the multi-pion production, and the 7r + production by similar contributions of 
all processes including the direct production. In this case, the A-resonance approximation is more 
accurate in terms of the shapes of the spectra. However, other processes quantitatively contribute 
as well. For the BB benchmark (lower panels), the tt production is dominated by the multi-pion 
production for most of the energies. Only in the low energy region, the charged pion production 
is governed by the direct and the resonant processes, and the 7r° production by the resonances. In 
this case, ignoring the high energy processes would lead to clearly misleading results. 

As far as the comparison among the different polarities is concerned, see Fig. [6J For all spectra, 
the 7r + are always most abundantly produced, as predicted above, followed by ir° and then 7r _ . For 
the GRB benchmark at high energies, where the multi-pion production dominates, the spectra are 
closest to each other. At low energies, there are significantly less tt~ produced than the other two 
polarities. However, note that, thanks to the multi-pion production, the tt~ are only suppressed 
by a factor of a few. The kaons, on the other hand, contribute about one to two orders less to 
the total meson fluxes. Nevertheless, there can be interesting effects in the high energy regime if 
cooling effects are present. For the AGN benchmark, because of the lower maximal proton energy 
times photon energy, even for large energies the tt~ are strongly suppressed. Otherwise, the result 
is qualitatively similar. For the BB spectrum we can nicely see the effect of lower multiplicities of 
7r~ compared to ir + and n° for high interaction energies (see Table [6|) in the spectrum at energies 
higher than 10 9 GeV. Otherwise, we have similar results as for the GRB benchmark. 



5. Comparison with SOPHIA 

Here we compare the results of our simplified models with each other and with SOPHIA. 
First, we focus on the primaries produced in the photohadronic interactions, mostly the pions. 
Then we compute and compare the neutrino spectra. In all cases, we use initial protons for the 
photohadronic interactions. 
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Figure 5 Contributions of resonant (thin solid), direct (dotted) and multi-pion (dashed) production for ty + , tt~ 
and 7T° spectra (from left to right) for proton-photon interactions. The upper panels are for the GRB benchmark, 
the middle ones for the AGN benchmark and the lower ones for the BB benchmark (see App. [C|. Computed with 
model Sim-B. 



-28 - 



GRB 



AGN 



BB 




10 2 10 3 10 4 10 5 10 6 10 7 10 8 10" 
£„/GeV 




10 3 10 4 10 5 10" 10 7 10 s 10" 10 10 
EJGeV 




10 4 10 5 10" 10 7 10 s 10 9 10 10 10 11 10 1: 
EJGeV 



Figure 6 Comparison among the 7r + (upper curve), tv° (middle curve), and 7r~ (lower curve) spectra for GRB 
(left), AGN (middle) and BB (right) benchmark. The grey curve shows in addition the K + spectrum. Computed 
with model Sim-B. 



5.1. Pion spectra 

The models considered in this subsection are listed in Table [7J where the (computational) com- 
plexity decreases from the top to the bottom. "SOPHIA" represents the output of the SOPHIA 
software, computed for our benchmark spectra from App. O The computation with SOPHIA is 
described in App. |Dl The model "BW" corresponds to the basic physics of SOPHIA as described 
in Sec. El including double integration for the direct production and multi-pion production from 
Sec. 14.4.21 In this case, we use Eqs. (|7j) and ([9]) for the secondary production using the two di- 
mensional response function including all resonances with full Breit-Wigner forms and the direct 
production as described in App. |A} Compared to SOPHIA, the kinematics is considerably simpli- 
fied. The multi-pion production is taken from Sec. 14.4.21 close to SOPHIA. "Sim- A" and "Sim-B" 
are the simplified models from Sec. using the factorized response function introduced in Sec. I4.H 
Whereas Sim-A treats all interaction types with the resonances explicitely, Sim-B defines an inter- 
action type for the A(1232)-resonance, and one for the higher resonances. Therefore, Sim-B uses 
considerably less interaction types. Note that we have obtained Sim-A and Sim-B by condensing 
the information from BW stepwise. Both models use the multi-pion production from Sec. 14.4.21 
whereas Sim-C is a simplified version of Sim-B including the multi-pion production in Sec. 14.4.11 

We compare the pion spectra produced by these different models in Fig. [7J (GRB benchmark), 
Fig- El (AGN benchmark) and Fig. [9] (Black body (BB) benchmark). In these figures, the upper rows 
show the pion spectra for 7r + , tt~ , and 7r° explicitly, whereas the lower rows show the pion ratios 
7r + /7r~, 7r + /7r°, and ir~ /it . Since the pion spectra for model BW, Sim-A and Sim-B are so close to 
each other that they can not be distinguished, we plot only the results of SOPHIA, model Sim-B 
and model Sim-C in the upper rows of Figs. [7] to El Whereas the charged pions lead to neutrino 
production, the neutral pions lead to photons. Therefore, the ratio 7r + /7r° determines, to leading 
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Abbrcv. Description 



Complexity 



Sec./Refs. 



SOPHIA SOPHIA software with full kinematics 

BW Resonances with full Breit-Wigner descrip- 

tion; direct production with full response 
function; simple kinematics 

Sim-A Simplified model with factorized response 
function and resonance treatment A (5- 
function approximation) 

Sim-B Simplified model with factorized response 
function and resonance treatment B (step 
function approximation) 

Sim-C Simplified model with factorized response 
function and resonance treatment B; sim- 
plified multi-pion production 



Monte Carlo method, 
~ 3000 xT Sim _ c 
Double integration, 45 IT, 
~ 120xT Sim _ c 

Single integration, 49 IT, 
~ (3-4)xT Sim _G 

Single integration, 23 IT, 
~ (2-3)xT Sim _ c 

Single integration, 10 IT, 

^Sim-C 



Miic ke et al 
(2000h 



Sec. [3] inch 
multi-pion 
from Sec. igj 
Sec. |U with 
Sec. EXT] and 
Sec. 14X21 
Sec. g] with 
Sec. 14X21 and 
Sec. 14X21 
Sec. H with 
Sec. 14X21 and 
Sec. HXTl 



Table 7 Considered models for the comparison of approaches. The (computational) complexity decreases from the 
top to the bottom. The time needed for the computation of the photohadronics in model Sim-C is given by Tsi m -c- 
The comparison of the computation times is done for power law spectra. For the computation with SOPHIA we used 
100000 trials per proton bin. 

order, the ratio between neutrinos and photons, which also often enters the computation of neutrino 
flux limits. On the other hand, the ratio tt + /ir~ affects the (electron) neutrino-antineutrino ratio. 
The ratio tt~ /it is shown for completeness. Note that the normalization of the different spectra is 
not chosen arbitrarily, but consistently to be able to compare the spectra directly; c/., App. iDl 

As the most important fact, as it can be read off from the upper rows of the three figures, all 
the spectra of our simplified models (apart from Sim-C, maybe) match the output from SOPHIA 
very well, both normalization and spectral shape. At high energies, however, where we imposed 
a sharp spectral cutoff, the spectra from SOPHIA are smeared out because of the more refined 
kinematics treatment, which can best be seen in Fig. [9] for the BB benchmark where we have 
a sharp spectral cutoff in the proton spectrum. This difference is unavoidable, and the price to 
pay for an efficient simplified model. Nevertheless, note that in more realistic spectra, or spectra 
averaged over different sources, such sharp features in the spectra may not be present. Although it 
seems that there are hardly any differences between SOPHIA and our models at lower energies, one 
can read off from the pion ratio plots in the lower rows (on a linear vertical scale) that there are 
small deviations. Whereas the differences at very low and high energies, where the spectra rapidly 
break off, are not very surprising, there are some differences coming from the different kinematics 
treatment. For example, SOPHIA actually produces even about 20% more ir + than -/r in the lower 
energy range for the GRB benchmark and the middle energy range for the AGN benchmark (c/., 
middle lower panels). We have checked that this difference neither comes from the resonance or 
direct production treatment, nor from the relative contribution of both processes. Instead, for some 
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Figure 7 Comparison of pion spectra (upper panel) and pion ratios (lower panel) for the different models from 
Table [7] for the GRB benchmark (see App.[Cj. Since the pion spectra for model BW, Sim-A and Sim-B can not be 
distinguished, we plot only the results of SOPHIA, model Sim-B and model Sim-C in the upper row. 



processes, our kinematics is a bit over-simplified, since for these a considerable amount of pions 
is (in SOPHIA) reconstructed at lower energies. Another example for an obvious difference is the 
high energy it + /ir~ and 7t + /tt difference in the lower panel of Fig. [9] for the BB benchmark, the 
most challenging one for the treatment of multi-pion production. The discrepancy is a result of the 
simplified kinematics treatment of taking the same x~ values for 7r + , 7r° and ir~ . This mainly effects 
the 7T~ spectra because they have, in comparison to ir + and 7r°, slightly different kinematics in 
SOPHIA. In the upper panel in Fig. [9] a double hump structure can be seen which follows from the 
kinematics treatment of multi-pion production that one part of the pions is reconstructed at lower 
energies (small x) an d the other at higher energies (large x)- If one averages over larger energy 
scales (such as in diffuse fluxes), such kinematics effects average out. 

As far as the comparison among our simplified models is concerned, the differences are small 
compared to the effects of kinematics discussed above. In fact, model "BW", which was originally 
designed as the most accurate reproduction of SOPHIA, produces the results farthest off from 
SOPHIA, especially for the lower half of the energy range. The reason may be that the errors 
introduced by the approximations in Sim-A and Sim-B partly compensate the errors from the 
simplified multi-pion production and direct channels. The model Sim-A was obtained from BW by 
assuming that all resonances are strongly peaked, whereas Sim-B was derived from this assumption 
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Figure 8 Comparison of pion spectra (upper panel) and ratios (lower panel) for the different models from Table[7]for 
the AGN benchmark (see App.[C|). Since the pion spectra for model BW, Sim-A and Sim-B can not be distinguished, 
we plot only the results of SOPHIA, model Sim-B and model Sim-C in the upper row. 



by collecting the properties of the higher resonances into one interaction type. In the comparison 
of model Sim-B to Sim-C the differences between the kinematics for multi-pion production from 
Sec. 14.4.21 and the simplified one from Sec. 14.4. II can nicely be seen. For the GRB (see Fig. [7]) and 
the AGN benchmark (see Fig. [8]) this mainly effects the high energy region whereas in the BB case 
(see Fig. [9]) most of the energy range is affected. Since the computation time for Sim-B, for which 
the results do not differ significantly from model Sim-A, is only about a factor of 2-3 longer and 
the high energy treatment is way more accurate than for Sim-C (especially close to the peaks), we 
focus in the following on model Sim-B. As one can see in Figs. [7] to El Sim-B is most accurate for 
power laws which are our main interest. Compared to SOPHIA, we gain a factor of about 1000 
(if implemented in C, depending on the integration method) in computation time for 100000 trials 
per proton bin in SOPHIA (as we use for the GRB benchmark). The Sim-B spectra do not have 
any small wiggles, because the computation is exact. However, note that the complexity of Sim-B 
increases with the number of interaction types, whereas the complexity of SOPHIA increases with 
the number of trials (and required smoothness of the functions). 
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Figure 9 Comparison of pion spectra (upper panel) and ratios (lower panel) for the different models from Table 
for the BB benchmark (see App.[C|. Since the pion spectra for model BW, Sim-A and Sim-B can not be distinguished, 
we plot only the results of SOPHIA, model Sim-B and model Sim-C in the upper row. 

5.2. Neutrino spectra 

In this section, we first review the production of neutrinos from pion and subsequent muon 
decays, as well as kaon decays. For the sake of completeness, we include the neutrinos from neutron 
decays, where the neutrons are produced by the photohadronic interactions. Then we compare our 
results to SOPHIA. We focus on model Sim B from the previous section only. 

The neutrinos are mostly produced in the following two decay chains: 



7T 



7T 



^ + -> e + + v e + , 
->■ (f~ + P M , 

\i~ -> e~ + v e + . 



(44) 
(45) 



For the energy spectrum from weak decays, we follow iLipari et al.1 (120071 ) (Sec. IV). In this case, the 
decays are described by Eq. ([2]). The functions dn i ^ b /dEi ) simplify, in a frame where the parent a 
is ultra-relativistic, to dn^ b / 'dEf, = 1/E a F a ^.b (Eb/E a ). T he functions F n _>h include the measured 
branching ratios in the possible final states and are given in lLipari et al.l (|2007l ) (Sec. IV). Note that 



-33- 



7r + and 7r are initially produced in different ratios and produce muons with different helicities, 
described by the scaling functions with r T = {rn^/m^) 2 : 

F ^$W = K-^-(-) = [l { ^)l Q{x ~ r - ] (46) 
F ^iV>= V-^*) = (T^^ 0(x - r - ) - (47) 

The muons decay further in a helicity-dependent way: 

F„ + ^{x,h)= F^-_ u ^x,-h) = f^-3x 2 + ^pj + /l (-^ + 3 ^-^r) ( 48 ) 
F^+^^x^h) = F^-^ e (x,-h) = (2 - 6x 2 + 4x 3 ) + h (2 - 12x + I8x 2 - 8x 3 ) (49) 

with h = 1 for right-handed and h = — 1 for left-handed muons. It is therefore mandatory to 
distinguish four muon states /j,^, [i£, and //^ as final states in order to account for the impact 
of muon polarization. The decay rates T?\ 6 = T in Eq. ([2]) (there is only one interaction type, 
which is decay) are just the inverse lifetimes T = r" 1 = Tq 1 m a E" 1 , where tq is the rest frame 
lifetime. 

For kaons, the leading decay mode into muon and neutrino is treated in the same way as in 



Lipari et al.l (120071 ) for the pion decays, i.e., with — >■ mx- The branching ratio for this channel 
is about 63.5%. The second-most-important decay mode is — > tt^ + 7r° (20.7%). The other 
decay modes account for 16%, no more than about 5% each. Since interesting effects can only be 
expected in the energy range with the most energetic neutrinos, we only use the direct decays from 
the leading mode. 

For protons accelerated in the jet, neutrons are produced by photohadronic interactions as 
described in App. |Bj Assuming that the neutrons escape from the acceleration region before they 
interact again, an additional neutrino flux from neutron decays is obtained. In this section, we show 
this neutrino flux separately. The beta decay describes the decay of the neutron into a proton, an 
electron and an electron anti-neutrino. In the ultra-relativist ic case, the mean fraction of the 



neutron energy going into the neutrino is x ~ 5.1 x 10 4 , see ( Lipari et al.l 120071 ). The neutrino 
spectrum is therefore obtained from the following equation: 

Qv e m = -Qn(—) (50) 
x V x J 

with Q n calculated from Eq. (|B6[) (App. [B]) . 

We show the neutrino spectra obtained from Sim-B and SOPHIA in Fig. [TUl where we also 
show the different contributions from different decay modes separately. Obviously, the SOPHIA 
and our combined spectra match very well, apart from the already discussed difference in the 
kinematics leading to some averaging in SOPHIA for large energies. Since the production of 7r + 
dominates for initial protons, Vn in Fig. [10] are most abundantly produced from pion decay; c/., 
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Figure 10 Comparison of spectra from the decays of different parents, as denoted in the labels. The left panel 
is for the GRB benchmark, the middle one for the AGN benchmark and the right one for the BB benchmark (see 
App.EJ. 




Figure 11 Comparison of the v £ (dashed) and v e (solid) spectra from the decays of different parents, as denoted 
in the labels. The left panel is for the GRB benchmark, the middle one for the AGN benchmark and the right one 
for the BB benchmark (see App. [C| . 

Eq. (|44p . However, the from muon decays, coming from the tt~ decay chain - cf., Eq. (|45p - are 
found at slightly higher energies and dominate for every high energies in the spectrum. For v^, the 
situation is exactly the opposite, but the final spectra look very similar. For very high energies, the 
SOPHIA spectrum is slightly higher than what one would expect, because other decay modes (such 
as from neutral kaons) contribute, which we have not considered. Without synchrotron cooling, 
the contribution from kaon decays is, however, small. 

In Fig. [TT1 we show the z/ e (dashed) and v e (solid) fluxes for our benchmarks. Obviously, the 
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Figure 12 Comparison of total electron to muon neutrino flavor ratio at the source for the following curves as 
given by the labels: Neutrinos from pion/muon decays, neutrinos from these including neutron and kaon decays, 
neutrinos from Sim-B (pion/muon decays) for without taking into account the spin state of the final muon (h — 0), 
and SOPHIA output. The horizontal lines mark the "standard" assumption for a flavor ratio electron to muon 
neutrinos 1:2. The left panel is for the GRB benchmark, the middle one for the AGN benchmark and the right one 
for the BB benchmark (see Add. [C| . 

v e , coming from the 7r + decay chain in Eq. (144h . dominate over the v e . However, if the neutrons 
produced by the photohadronic interactions escape from the source and then decay, they will lead 
to an additional neutrino flux shown by the thin gray curves (not included in the total v e curves). 
Especially at very low energies, the v e flux then dominates. 



5.3. Flavor and neutrino-antineutrino ratios of the neutrinos 

In this section, we discuss the electron to muon neutrino flavor ratio (the ratio between the 
electron and muon neutrino fluxes) and the neutrino-antineutrino ratios at the source. The flavor 
composition at the source is primarily characterized by the flux ratio (p e + i/ e )/(i/ /J + v^), since 
almost no v T (or v T ) are expected to be produced at the source. Because neutrino telescopes can, 
in muon tracks or showers, not distinguish neutrinos from antineutrinos, this ratio is representative 
for the detection as well. From Eqs. (j44]) and (|45|) . we can read off that this flavor ratio should 
be about 1/2 without kinematical effects, which we call the "standard assumption". At the Earth, 
the thr ee neutrino flavors then a lmost equally mix (in the ratio 1:1:1) through neutrino flavor 



mixing ([Learned fe: Pakvasalll995l ). 



We show in Fig.[l2]the flavor ratio at the source for the GRB (left panel), AGN (middle panel) 
and BB (right panel) benchmark. The standard assumption 1/2 is marked by the horizontal lines. 
Our curves from Sim-B (thick solid, neutrinos from pion and muon decays only) bend upward from 
this standard assumption, whereas the SOPHIA curves bend downward for large energies. This 
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Figure 13 Comparison of the neutrino-antineutrino ratios at the source for electron neutrinos (left panels) and 
muon neutrinos (right panels) for the following curves given by the labels: Neutrinos from pion/muon decays, neutrinos 
from these including neutron and kaon decays, and SOPHIA output. The horizontal lines in the lower panels mark 
the "standard" assumption for a flavor ratio muon neutrinos to antineutrinos 1:1. The left panels are for the GRB 
benchmark, the middle ones for the AGN benchmark and the right ones for the BB benchmark (see Add. [C]) . 

difference can be explained by a different implementation of the weak decays: if we do not take 
into account the spin state of the final muon (dotted curves h = 0), we can reproduce the SOPHIA 
results almost exactly with Sim-B. In fact, the effect of the helicity is larger than the details of the 
interaction model. The dashed curves include the effect of neutron decays (low energies) and kaon 
decays (high energies) into Sim-B. Especially for low energies, where v e are abundantly produced 
by neutron decays, the curves deviate from the thick ones. This effect is strongest for the AGN 
benchmark, for which the standard assumption 1 /2 only approximately holds in a relatively narrow 
energy window. Note that the dashed curve i n the left panel and all the other results for the GRB 
benchmark exactly match Lipari et al. (j2007i ). where the weak decays are discussed in detail. 



If the Glashow resonance process v e +e -^-W — > anything at around 6.3 PeV can be observed 
in a neutrino telescope, the neutrino-antineutrino ratios at the source may be relevant as well (all fla- 
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vors at the source contribute to the flux at the Earth throu gh flavor mixing) (jLearned <fe Pakvasa 



19951 ; lAnchordoqui et al.ll2005l ; lBhattachariee &: Guptall2005l ). The neutrino-antineutrino ratio may 



be relevant to distinguish p7 interactions at the source, for which mostly ir + are produced, from pp 
interactions at the source, for which ir + and ir~ are produced in almost equal ratios. Therefore, we 
show in Fig. [I3]the neutrino-antineutrino ratios at the source for electron neutrinos (left panels) and 
muon neutrinos (right panels). The electron neutrino-antineutrino ratios in the left panels depend 
on the ratio of ir + and ir~ produced, see Eqs. (|44p and (j45]). Our result matches SOPHIA very 
well, especially in the important energy range from the peak of the spectrum two decades down, 
apart from the discrepancy for high energies for the BB benchmark (upper right panel) which we 
discussed already as it is coming from the tt + /tt~ ratio. The deviation between SOPHIA and Sim-B 
can be up to 30%. After flavor mixing, the correction to the electron neutrino-antineutrino ratio 
at the Earth is at the level of 10%, much smaller than the effect on the flavor mix expected from 
pp interactions. The muon neutrino-antineutrino ratios in the right panels do not depend on the 
ratio of tt + and 7r~ produced, as in every pion decay the same number of muon neutrinos and 
antineutrinos is produced; see Eqs. (|44p and (|45p . In this case, our results match SOPHIA and the 
standard prediction very well, and the effects of neutron and kaon decays are small in the absence 
of cooling. 



6. Summary and conclusions 

We have discussed simplified models for photo-meson production in cosmic accelerators. The 
main purpose of this simplification has been the definition of a photohadronic interaction model 
useful for efficient modern time-dependent AGN and GRB simulations, and for large-scale param- 
eter studies, such as of neutrino flavor ratios. The major requirements have been listed in the 
introduction. For example, the secondaries (pions, kaons) are not to be integrated out, since their 
synchrotron cooling affects the neutrino flavor ratios. 

We have first re-phrased the problem in terms of a two-dimensional "response function" to 
be folded with arbitrary photon and proton input spectra in order to compute the secondary 
fluxes. The key idea for our simplified models has been the factorization of this two-dimensional 
response function, which has allowed to eliminate one of the integrations. In order to include 
kinematics as good as possible, we have then defined a discrete number of different interaction 
types with different characteristics based on the underlying physics of SOPHIA. The kinematics 
of more complicated interactions, such as direct production or multi-pion production, has been 
simulated by the introduction of multiple interaction types for each production channel. In a step- 
by-step fashion we have simplified then the resonance treatment in order to arrive at our simplified 
model Sim-B. It allows for the computation of pion spectra with only one integral, summed over 
about ten interaction types, and can be easily adopted from our description. The extendibility of 
this approach has been demonstrated by showing how K + fluxes can be added, once a suitable 
parameterization is found. 
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Similarly, the response function can be easily changed if new data are provided, new processes 
can be included, or systematics on the particle physics can be added. Of course, some effort has to 
be spent to find a suitable parameterization for each process. 

We have demonstrated that our results match the output of SOPHIA sufficiently well. How- 
ever, there are some differences due to the more refined kinematics treatment of SOPHIA, which 
effectively corresponds to one additional integration. For example, for very narrow spectral features, 
such as rapid cutoffs, the spectra are naturally more smeared out by SOPHIA. This is especially 
the case for high energetic interactions where multi-pion production is dominant, as can be seen 
in the BB (black body) benchmark. However, our approach is much simpler in the sense that the 
interaction rate and the folding with the proton spectra is automatically taken into account (c/., 
App.[Dj). In addition, we hay e included the spin state of the final muon in the pion decays, as de- 



scribed in lLipari et al.1 (|2007l ). not included in SOPHIA, which leads to differences in the neutrino 
flavor ratios: in fact, the electron to muon neutrino flavor ratio at the source is typically larger 
than 0.5, instead of smaller, as predicted without taking into account the spin state. In particular, 
we obtain v e : v„ : v T ~ 1 : 1.85 : for the GRB benchmark, 1 : 1.96 : for the AGN benchmark, 
and up to 1 : 1.82 : for the BB benchmark close to the spectral peaks. This means that especially 
the AGN benchmark behaves as pion beam in spite of the helicity dependence of the muon decays, 
whereas the BB benchmark shows the strongest deviation. 

Since our approach has allowed us to discuss the leading interaction types separately, w e have 



also shown the differences to the A(1232)-resonance approximation (see also iMiicke et al.1 (|2000l . 



19991 )). For example, we have shown how multi-pion production modifies the characteristic shape 
of the GRB neutrino spectrum expected from the resonance approximation (which is flat in E 2 
times the flux in the intermediate energy range). In fact, all of the the resonances combined do 
not dominate in any significant part of the charged pion spectrum for our GRB, AGN and BB 
benchmarks. In addition, the A(1232)-resonance approximation has rendered insufficient for the 
computation of the (electron) neutrino-antineutrino ratios at the source, because, by definition, 
no 7r~ are produced. We have also demonstrated from our general response function in model 
Sim-B, that for any input proton or photon spectrum the tr + /ir° ratio at the source is significantly 
larger than one, as opposed to 1/2 from the A(1232)-resonance approximation. This implies that 
any neutrino flux based on the A(1232)-resonance approximation and normalized to photon flux 
observations using the charged to neutral pion ratio is underestimated by at least a factor of 2.4, 
where this factor is independent of the input spectra. 

We conclude that our simplified model Sim-B allows for an efficient computation of pion and 
neutrino fluxes, including the necessary features for neutrino flux ratio discussion and the necessary 
efficiency for time-dependent simulations. It is sufficient for many purposes especially for power 
law photon fields, but, of course, it cannot replace a full Monte Carlo simulation including full 
kinematics if high precision fits of existing data are required. In particular, it may turn out to be 
useful for time-dependent simulations and extensive parameter space studies using power law spec- 
tra, including spectral breaks. However, in other cases, such as if the high energy interactions with 
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photon spectra with sharp peaks are very important, or there are anisotropic photon distributions, 
the Monte Carlo method may be the better choice. 
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A. Kinematics treatment of direct production 
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Figure 14 Energy fraction \ g°i n g into the pion as a function of e r for direct one pion production for different 
values of cos 6 and our simplified approach "SIM" ; see main text for details 



Here we follow the approach of lRachenl (|1996l ). As mentioned in Sec. 13.21 the direct production 
process is strongly backward peaked with respect to the produced pion. One possible approximation 
for the determination of x is to assume that the average scattering angle is 180°. The energy fraction 
X going into the pion as a function of e r is shown in Fig. [14] as the curve cos# = —1 (c/., Eq. (I16p ). 
For comparison, the curves cos# = and cos# = 1 are shown as well. For a more accurate 
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represen tation of kinematic s, we take the probability distribution of the Mandelstam variable t, as 
given in 



Miicke et al.l (|2000l ). into account. For small \t\, it is given by: 

~dt = b e bt (~ 1 ) -e bt W 



(Al) 



with b ~ 12GeV 2 . For example, for the interaction type Ti, the Mandelstam variable t is given 
by 

t{ C o^) = ml-2 S ^ S + m ^ ml - 2 -^\\[ - 1 7^ ) -m^cos^. (A2) 




Combining Eqs. (|A1|) and (|A2p . we obtain the average scattering angle as a function of the center- 
of-mass energy. Inserting the result into Eq. (|16p . we obtain the average energy fraction \ g° m g 
into the pion for direct one pion production. It is shown in Fig. [TH as curve "< cos 9 >" as a 
function of e r . Analogously we compute the mean % f° r direct two pion production by combining 
Eq. (jAip with the variable t for the considered process. 



In the simplified model for direct production (c/., Sec. 14. 3p we use the factorized response 
function (see Eq. (|27p ) as for the simplified models of the other processes. Therefore we have to 
choose an energy independent, constant x- Since the range of x-values for direct production is 
wide, we divide the e r -range into three sections (for each of the interaction types Ti and T2 a ) 5 such 
that it reproduces the results of the non-simplified model for typical power law spectra for photons 
and protons in astrophysical sources, such as GRBs and AGNs. This approach corresponds to the 
thick gray curve "SIM" in Fig. [T4"l where only the lower two interaction types are visible in the 
plotted energy range. Obviously, it is a good step function approximation to the curve "< cos 6 >" . 
The different sections of this step function correspond to our interaction types "L" , "M" , and "H" . 



B. Cooling and escape of primaries, re-injection 

A related issue to the secondary particle production is the cooling or escape of the primaries 
due to the interaction process. We do not focus on the cooling or escape timescales in this study, 
but, for the sake of completeness, we demonstrate how they can be computed from the quantities 
presented for our simplified models in Sec. [H If the primaries lose energy in an interaction, such 
as protons or neutrons in pion photoproduction p + 7 — > p + 7r°, this process can be interpreted as 
cooling, whereas if the primaries disappear, such as protons in p + j — > n + 7r + , it can be interpreted 
as escape. In the latter case, neutrons are re-injected into the system, which we will discuss below. 
The cooling rate t~ ool (E p ) = —E~ 1 dE p /dt or escape rate = —N~ 1 dN p /dt for the species p 
(proton or neutron) due to the photohadronic interactions is, for constant -fT IT ', given by 

d(^) = E r P T >p(^) ifIT > E f1 ^A e p)- ( B1 ) 

IT it.pVp 
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in terms of the quantities in Eq. ([2]). Here K 1 ^ ■ E p is the loss of energy per interaction; therefore, 
if IT is often called "inelasticity". Note that if it is a function of the kinematical variables, such as 
the center of mass energy, it has to be folded into the calculation of the interaction rate in t~\. 
However, in Sec. [U we have constructed the interaction types such that X IT is a constant. For 
photo-pion production, the inelasticity can be related to the x l J-^b ™ Eq- © ^ 

K IT = J2^ b Mf, (B2) 

i.e., the energy loss of the nucleon equals the energy deposited in all interaction products (other than 
the initial nucleon). Note that the classification as cooling or escape also depends on if protons and 
neutrons are distinguished in the final state. In this section, we distinguish protons and neutrons. 
In addition, note that in astrophysical objects there may be other sources of cooling and escape to 
be taken into account, such as synchrotron cooling or escape from the production region. 

The quantity needed for the computation of Eq. (|Bip is the interaction rate in Eq. Q . Compar- 
ing Eq. (|4|) with Eq. (127p . we find for our simplified models that the interaction rate for interaction 
type IT of the initial nucleon p is 

oo 

T IT (£ P )= J dsn^f^^f) , (B3) 



i.e., conveniently parameterized in terms of our simplified response function. Then the cooling 
and escape rates in Eq. (|Bip can be written in terms of the initial M p or different nucleon M p i 
multiplicity (M p + M p , = 1): 

Ci(^) = E< rIT (^)^ IT > *«W= E <r IT (£ p) . (B4) 

IT IT, pVp 

The nucleon multiplicities are for the resonances in model A given in Table [T] for the interaction 
types in Table [3l for the resonances in model B in Table HI for multi-pion production M p=p i = 0.69 
and M p ^ p i = 0.31, and for direct production in Table [1] for the interaction types in Table [5j The 
inelasticities can in resonance model A be obtained from the x IT i n Table [3] according to Eq. (|B2j) , 
i.e., K 1T 1,2,3 = J2n Xp^-vr -^ T > where the number of pions produced M^ T = 1 for ITs 1, 2a, and 2b, 
and Ml T = 2 for IT 3 (for IT 2, the pions in 2a and 2b are summed over). For resonance model B, 
they are given in Table HI for multi-pion production x Mnltl ~ n ~ 0.6, and for direct production they 
are listed in Table [5] (here IT T2b is not counted separately) . 

The re-injection rate p — > p 1 for initial nucleons p can be obtained analogously to Eq. (|27h and 
Eq. ([25} from 

R lT (x, y) = 5{x - (1 - K IT )) M l p 7 f T (y) (B5) 



as 



oo 
eth/2 
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Note that double counting of the same interaction has to be avoided. In particular, interaction 
type 2 must not be counted twiceQ In addition, note that for multi-pion production Eq. ([38]) based 
on the total cross section should be used in all cooling, escape, and re-injection rates. 



C. Benchmarks 



Our benchmarks are given in the SRF. The benchmark for GRBs taken from iLipari et al 



(|2007l ) . The photon spectrum, a broken power law, is given by 

f N 7 (gIv)" 1 0.2eV< 



n 7 (e) = 
and the proton spectrum by 



_N 7 10" 6 



(cev) 



GeV cm^ 

2 1_ 

GeV cm 3 



IkeV < e < 300 keV 



(CI) 



n p {E p ) = N p (tJ^J exp 



6.9 • 10 8 GeV 



1 



GeV cm 3 



Ep > 1 GeV. (C2) 



Note that there are dimensionless normalization constants N 7 and N p . The resulting neutrino 
spectrum is characterized by a wide maximum in E 2 Q U (E). This benchmark is designed to fit 
an average of the total distribution of GRB. We limit ourselves to this average, since taking all 
extreme cases of GRB would go beyond the scope of this paper. 



The benchmark for AGNs is adopted from iMiicke Sz Protherod (120011 ). The photon spectrum, 
a power law with a sharp cutoff, is given by 



n 7 (e) 



'n ( £ 1 

1N 7 \GeVJ GeV cm 3 



:i-4-io- 7 ) - 2 (cev)" 1 ' 



1 

GeV cm 3 



10~ 3 eV < e < 140 eV 
140 eV < e < 3.6 keV 



and the proton spectrum is given by 

n p {E p ) = N p (-^) exp 



E n 



2.6 • 10 9 m t 



1 



GeV cm 3 



E p > 1 GeV 



(C3) 



(C4) 



with m p = 0.938 GeV. This benchmark is well in the range of usual parameters of the HBL 
(a subclass of AGN), which are the most interesting objects for state-of-the art Air Cerenkov 
Telescopes. 

The thir d benchmark, t he mo st challenging one, is a high energetic black body spectrum (BB), 
adopted from lMiicke et al.l (|1999l ). The photon spectrum of temperature 10 eV is given by 

nJe) = N 7 1.318 • 10 31 (7^77) 2 r— ^ 7 

7W 7 VGeV/ exp • 10 8 ] - 1 GeV cm 3 



(C5) 



7 Although there are two pions produced, there is only one secondary nucleon. For the inelasticity, however, the 
energy losses into all pions have to be taken into account. For the interaction rate, ITs 2a and 2b are counted as one. 
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and the proton spectrum with a sharp cutoff is 
n p (E p 



iN P I GcV 



GeV cm 







10 6 GeV < E p < 10 12 GeV 
else. 



(C6) 



The blackbody temperature is designed to fit usual BLR photons. It is the most challenging bench- 
mark because the proton spectrum is not smeared out for the highest energies by an exponential 
cutoff and the photohadronic interactions are dominated by high energetic multi-pion production 
due to the conside red peaked photon spectrum. Even though the benchmark is adopted from 
Miicke et al.l (|1999l ). it does not exactly represent real physics. Thermal spectra are usual produced 
outside the shock and are therefore b eamed. Only for the production of cosmogenic neutrinos with 
CMB photons beaming is negligible (|Steckerlll979l ) 



D. Comparison with SOPHIA runs 

In SOPHIA, an injected proton of energy E p is assumed to interact for sure, and the secondary 
particle (of type b) distribution dn p ^.b/dEb(E p , Eb) is computed for a specific proton energy or a 
range of proton energies. This means that 



/ 



^L(E p ,Eb)dEb = N b (Dl) 
aEb 



is the number of secondary particles produced. The particle spectrum can then be computed using 

Qb(E b ) = J dE p Np(Ep) T p ^ b (E p ) ^j^(E p , E b ) . (D2) 
Note that this formula is very similar to Eq. ([2]), but not split up in different interaction types. 



As it is obvious from Eq. (jD2[) , the interaction rate as a function of proton energy is needed as 
additional input. We use Eq. with the total cross section as depicted in Fig. Q] and parameterized 
in Sec. 13.31 to compute the interaction rate (the cross section is already summed over all interaction 
types then). As far as the units and normalization are concerned, it is useful to specify all energies 
in GeV and all cross sections in /ibarn = 10 -30 cm 2 (note, however, that in SOPHIA, photon spectra 
are always given in eV). In this case, the interaction rate carries units of 3 • 1(T 20 N 7 s~\ where N 7 
is the dimensionless normalization of the photon spectrum in App.O From Eq. (|D2|) . it is obvious 
that Qb comes in units of 3 • 10 -20 N 7 N p s _1 GeV -1 cm~ 3 if N p is given in units of GeV -1 cm -3 
and carries the (dimensionless) normalization factor N p . The same units apply to the results from 
our simplified models. 

SOPHIA uses logarithmic energy spacing in E p and Eb and provides the output on a discretized 
energy grid in the form dn p ^b/dxb(E p , Xb) with Xb = \og{Eb/E p ), equally spaced in Ax p with 
x p = log(E p /GeV) and in Axb- For the easiest data extraction, it is advisable to use Ax p = Axb- 
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Then, it is useful to re-write Eq. (|D2|) as 

Q b ( x {) = A Xp n p (4) r ^K) ^^(4> 4) w ~ xi . (° 3 ) 

where the last factor comes from switching to the logarithmic scale. With Eq. (|D3j) . the SOPHIA 
output can be used directly, collecting all entries for a specific x 3 h from all proton energy bins. Note, 
however, that still the particular output format has to be taken into account (SOPHIA first lists 
the bin range filled with data, then the filled bins, as a function of the proton energy). 

For the test runs of the AGN and GRB spectra, we use a proton energy grid between 1 GeV 
and 10 10 GeV with 100 bins, i.e., Ax p = 0.1. In addition, we use 100 output bins with a step size 
Axb = 0.1. For the GRB benchmark, we use 100000 trials per proton bin, for the AGN benchmark 
25000 trials. For the BB test run, we use a proton energy grid between 10 6 GeV and 10 12 GeV with 
60 bins, i.e., Ax p = 0.1. In addition, we use 75 output bins with a step size Ax b = 0.1 and 10000 
trials per proton bin. 
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